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We have developed a set of techniques for performing large 
scale ab initio calculations using multigrid accelerations and a 
real-space grid as a basis. The multigrid methods permit effi- 
cient calculations on ill-conditioned systems with long length 
scales or high energy cutoffs. The technique has been applied 
to systems containing up to 100 atoms, including a highly 
elongated diamond cell, an isolated Ceo molecule, and a 32- 
atom cell of GaN with the Ga d-states in valence. The method 
is well suited for implementation on both vector and massively 
parallel architectures. 



Over the course of the last several decades, plane- 
wave-based methods have been used to perform elec- 
tronic structure calculations on a wide range of phys- 
ical systems. The Car-Parrinello (CP) and other it- 
erative methodalrEl have made such calculations pos= 
sible for systems containing several hundred atoms.u 
While these methods have been very successful, sev- 
eral difficulties arise when they are extended to phys- 
ical systems with large length scales or containing 
first-row or transition-metal atoms. Special techniques 
have been developed to handle soirie of these prob- 
lems. Optimiaed pseudopotentialspQ the augmented-, 
wave method,0 plane waves in adaptive-coordinates,cl 
and preco nditi oning combined with conjugate-gradient 
techniqueal3'llll have had considerable success. How- 
ever, these techniques are still constrained by the plane- 
wave basis set, which requires periodic boundary con- 
ditions for every system and fast Fourier transforms 
(FFTs) to efficiently transform between real and recipro- 
cal space. FFTs involve non-local operations that impose 
constraints on the adaptability of these algorithms to 
massively parallel computer architectures, because they 
perform best on problems that can be divided into local- 
ized domains. 

It has been appreciated for some time that there are 
potential advantages to performing electronic-structure 
calculations entirely in real space. Boundary conditions 
are not constrained to be periodic, which permits the use 
of non-periodic boundary conditions for clusters and a 
combination of periodic and non-periodic boundary con- 
ditions for surfaces. By employing nonuniform real-space 
grids, it is possible to add resolution locally; e.g., for a 
surface or cluster calculation, a basis that uses a high 
density of grid points in the regions where the ions are 
located and a lower density of points in the vacuum re- 
gions can lead to order of magnitude savings in the basis 



size and total computational effort £3 More importantly, 
the use of a real space basis opens up the possibility of 
using multigrid iterative techniques to obtain solutions of 
the Kohn-Sham equations. Multigrid methodal3 provide 
automatic preconditioning on all length scales, which can 
greatly reduce the number of iterations needed to con- 
verge the electronic wavefunctions. Furthermore, the 
real-space multigrid formulation does not involve long- 
range operations and is particularly suitable for pax=. 
allelization and wavefunction-based 0(N) algorithms, Ed 
because every operation can be partitioned into hierar- 
chical real-space domains. In this Communication we 
describe the essential elements of the uniform-grid, real- 
space multigrid implementation as well as tests on a 
variety of systems, including a vacancy in diamond, a 
strongly elongated 96-atom diamond supercell, and a 32- 
atom supercell of GaN with 3d electrons in valence. 

There have been a number of previous real-space 
grid-based electronic-structure calculations. Tie finite- 
element method was applied by White et alx3 to one- 
electron systems. They used both conjugate-gradient 
and multigrid acceleration to find the_ground-state wave- 
function. Two of the present authorsEj used nonuniform 
grids with locally enhanced regions in conjunction with 
multigrid acceleration to calculate the electronic proper- 
ties of atomic and diatomic systems with nearly singular 
pseudopotentials. They verified that the preconditioning 
afforded by multigrid was effective injnulti-length-scale 
systems. Recently, Chelikowsky et al&Q have used high- 
order finite-difference methods and soft non-local pseu- 
dopotentials on uniform grids to calculate the electronic 
structure, geometry, and short-time dynamics of small Si 
clusters. 

Several issues that are absent from plane-wave or 
orbital-based methods arise when using a real-space grid 
approach. In the former case the wavefunctions, poten- 
tials, and the electronic density are representable in ex- 
plicit basis functions, and thus are known everywhere. 
Errors in the representation of these quantities are mainly 
due to the truncation of the basis. In a real-space grid 
implementation, these quantities are known only at a dis- 
crete set of grid points, which can introduce a spurious 
dependence of the Kohn-Sham eigenvalues, the total en- 
ergy, and the ionic forces on the positions of the ions 
with respect to the real-space grid. We have developed a 
set of techniques that can overcome these difficulties and 
can be used to compute accurate static and dynamical 
properties of large physical systems, while taking advan- 
tage of the rapid convergence rates that are made pos- 
sible by multigrid methods. In our formalism the wave- 
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functions, density, and potentials are directly represented 
on a uniform three-dimensional real-space grid with lin- 
ear spacing h gr id and number of mesh points Ngrid- 
The ions are represented by soft-core norm-conserving 
pseudopotentialsO Exchange and correlation effects are 
treated using the local density approximation (LDA) of 
density functional theory. 

We start with the eigenvalue expression for the total 
electronic energy of a system of electrons and ions in the 
local density approximation (LDA): 

Eld A = ^2 6,1 + / d ?p{r) (e X c(p(r)) ~ Pxc(p(r))) 

\ J dfr P^ VHartree ^ + E lon ^ ion . (1) 

Since the kinetic energy operator is not explicitly rep- 
resented in Eq. ([!]), we can discretize the Kohn-Sham 
equations in a generalized eigenvalue formulation: 

A[i/) n ] + B[V eff tp n ] = e n B[ip n }. (2) 

A and B are. the components of the Mehrstellen 
discretization,II3 which is based on Hermite's generaliza- 
tion of Taylor's theorem. It uses a weighted sum of the 
wavefunction and potential values to improve accuracy of 
the discretization of the entire differential equation, not 
just the kinetic energy operator. The weights are listed 
in Table ||. The Kohn-Sham hamiltonian, wavefunctions, 
and eigenvalues are represented to 0(h grid ), but the pref- 
actor for the 0(hg rid ) term is smaller than 0.01. Alterna- 
tive grid-based discretizations of the Kohn-Sham equa- 
tions do exist, e.g., high-order central-finite-difference 
methods were used by Chelikowsky et alX3 They dis- 
cretize the kinetic energy operator only, and typically use 
an 0(h g ^. id ) discretization. Mehrstellen discretization dif- 
fers from central-finite-difference discretization in that it 
achieves higher accuracy by using more local data. For 
example, the fourth-order A operator extends to second 
nearest neighbors; in contrast, the twelfth-order central- 
finite-difference operator extends linearly out 6 mesh 
points in each direction. We have found in our test cases 
that the fourth-order Mehrstellen discretization produces 
equivalent or better accuracy than the standard sixth- 
order finite-difference discretization. 

When using plane- wave basis sets, the accuracy is usu- 
ally determined by the convergence of the total energies 
and eigenvalues as the energy cutoff used is increased. In 
calculations employing a real-space grid, it is also nec- 
essary to consider the dependence of the energies and 
the ionic forces on the positions of the ions relative to 
the real-space grid points. This non-physical dependence 
is not present in plane- wave methods. Our selection of 
the Mehrstellen operator and discretization techniques is 
designed to reduce this grid dependence as much as pos- 
sible, because the determination of accurate trajectories 
in molecular dynamics simulations and the calculation 



of dynamical quantities is very sensitive to this depen- 
dence. Our approach is to perform calculations on iso- 
lated atoms and to monitor the variations as the ion is 
moved relative to the grid points. To reduce finite-size 
effects, periodic boundary conditions and a large super- 
cell are used in these tests rather than cluster boundary 
conditions. In Fig. [I] we show the variation in the total 
energy of a carbon atom. For bulk calculations, we also 
test the grid resolution by monitoring the total energy 
as the atomic lattice is rigidly translated with respect to 
the grid. 

The solution of Eq. (^) by direct matrix methods or 
by standard iterative techniques is prohibitively expen- 
sive for very large systems because they require 0(N^f d ) 
operations per wavefunction on an uniform mesh of N gr id 
pointsJ13 As the grid resolution is increased, e.g., to treat 
systems containing first row or transition metal atoms, 
the rate of convergence becomes unacceptable. A sim- 
ilar slow-down is present in plane-wave methods as the 
energy resolution increases. 

To efficiently solve Eq. (0), we have used multigrid 
iteration techniques that accelerate convergence by em- 
ploying a sequence of grids of varying resolutions. The 
final solution is obtained on a grid fine enough to accu- 
rately represent the pseudopotentials and the electronic 
wavefunctions. If the solution error is expanded in a 
Fourier series, it may be shown that iterations on any 
given grid level will quickly reduce the components of 
the error with wavelengths comparable to the grid spac- 
ing but are ineffective in reducing the components, with 
wavelengths large relative to the grid spacing.E3o The 
solution is to treat the lower frequency components on a 
sequence of auxiliary grids with progressively larger grid 
spacings, where errors appear as high frequency compo- 
nents. This procedure provides excellent preconditioning 
for all length scales present in a system and leads to very 
rapid convergence rates. The operation count to converge 
one wavefunction with a fixed potential is 0(N gr id), com^ 
pared to 0(N gri d logN gri d) for FFT-based approaches.H 
In addition, all operations except orthogonalization are 
short ranged, which allows for easy parallelization and a 
natural implementation of O(N) algorithms l3 

One iteration towards the solution of Eq. (|l|) consists 
of a multigrid step to solve Eq. (Q), followed by orthog- 
onalization of the orbitals, and an update of the elec- 
tronic density. The components of the effective poten- 
tial that depend on the density (Hartree and exchange- 
correlation) are then recomputed for the new density and 
the process is repeated until self-consistency is reached. 
We have found that it is often necessary to implement a 
mixing procedure when updating the electronic density to 
minimize charge sloshing .o The Hartree potential is com- 
puted by solving Poisson's equation using multigrid iter- 
ations on the corresponding Mehrstellen discretization. 
The non-local ionic pseudopotenti&Uis evaluated in real 
space using the Kleinman-BylanderEil separable form. To 
facilitate comparisons with plane- wave calculations, we 
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define an equivalent energy cutoff for the multigrid cal- 
culation, 7r 2 /2h 2 [Ry], to be equal to that of a plane- wave 
calculation that uses a FFT-ffrid with the same spacing 
as the multigrid calculation^ Using this convention, the 
computational time required to perform one step of the 
above procedure is approximately the same as for a single 
step in the CP method when the cutoffs are equal. 

The multigrid method outlined above has been tested 
on several large systems. A 64-atom diamond super- 
cell was chosen as a representative periodic system for 
which both multigrid and Car-Parrincllo calculations are 
feasible. The calculations were performed on the per- 
fect crystal and— the relaxed vacancy using the same 
pseudopotential.O The CP cutoff was 35 Ry, while the 
multigrid code used a grid spacing of 0.336 bohr, which 
results in a grid with the same total number of points 
as the FFT grid needed in the CP code. The k-space 
sampling was restricted to the T point. The dynamical 
relaxation methooa was used to quickly relax the va- 
cancy in both calculations. 

The Car-Parrinello and real-space calculations are in 
excellent agreement with each other. The relaxed ionic- 
positions agree to within 0.009 bohr for all ions in the 
supercell, and the largest difference in Kohn-Sham eigen- 
value is 0.06 eV. The diamond vacancy introduces an 
initially triply degenerate level in the gap, which splits 
into a singlet and a doublet due to the Jahn- Teller ef- 
fect. The magnitude of this splitting is 0.32 eV in both 
the real-space and Car-Parrinello results. The main re- 
sults are summarized in Table [H]. The discrepancies be- 
tween the LDA calculations and experimental values for 
the bandgap and the cohesive energy are well-known de- 
ficiencies of density functional theory. 

An isolated Ceo molecule was selected as an example 
of a non-periodic system. The simulation cell was a cube 
of length 23 bohr and the grid spacing was 0.360 bohr. 
The initial ionic coordinates were.-tfenerated using the 
classical Tersoff-Brenner potential,E-3 and the electronic 
wavefunctions were set to zero on the boundaries of the 
cell. After the convergence of the electronic system, the 
ions were relaxed using the same relaxation scheme as be- 
fore. Two distinct bond lengths were found in the final 
structure, corresponding to the carbon-carbon single and 
double bonds. There were twice as many single bonds as 
double bonds, and the average double and single bond 
lengths were 1.39 and 1.44 A, compared to-4.41 and 1.45 
A obtained in a previous CP calculations for the Ceo 
solid. The standard deviations of the bond-length distri- 
butions were on the order of 10~ 3 A in both calculations. 
The experimental values for the solid are 1.40 and 1.45 
A, respectively. 

A significant advantage of the multigrid method is the 
speed of convergence to the electronic ground state for 
a given initial ionic configuration. For systems requiring 
a small energy cutoff or of small size, the speed advan- 
tage with respect to CP-based methods is not substan- 
tial. However, for systems requiring a large energy cutoff, 
or of large dimensions, this advantage — as measured in 



actual computational time — is typically an order of mag- 
nitude. This is because the maximum stable timestep in 
the Car-Parrinello method must be much smaller than in 
the multigrid approach. 

To illustrate the ability of the multigrid method to 
handle ill-conditioned systems, test calculations were per- 
formed on a 32-atom GaN cell in the zinc-blende struc- 
ture with the Ga d-electrons in valence, and on a highly 
elongated 96-atom diamond cell of dimensions (6.72, 
6.72, 80.64) bohr. For the GaN cell a grid spacing of 0.175 
bohr was used, which corresponds to an energy cutoff of 
160 Ry in a plane-wave calculation. Only the T point 
was included in k-space sampling. Starting with random 
initial wavefunctions, 240 self-consistent steps were re- 
quired to converge the total energy. Recently, several 
calculations have been performed on GaN that explicitly 
included the d-electrons in valence; the multigrid results 
are in good agreement with these calculations (see Ta- 
ble III). The calculations for the elongated diamond cell 
used a grid resolution of 0.336 bohr, corresponding to an 
energy cutoff of 35 Ry in a plane wave calculation. For 
convenience, an approximate solution was generated on a 
grid with twice the final spacing, on which each iteration 
is eight times faster than on the final grid. After transfer- 
ring this approximate solution to the final grid, only 50 
self-consistent iterations were sufficient to converge the 
total energy to a relative tolerance of 10~ 8 . 

These two ill-conditioned systems represent worst-case 
scenarios, and the performance of the multigrid method 
is considerably better for typical systems. Table |v| il- 
lustrates the convergence properties as a function of grid 
resolution for an 8-atom diamond cell. The observed con- 
vergence rates are largely independent of the energy cut- 
off. The number of self-consistent steps required to con- 
verge the density is also nearly independent of the sys- 
tem size. At an equivalent cutoff of 35 Ry, the multigrid 
method required 17 steps to converge the total energy of 
the 8-atom diamond cell to a tolerance of 10 -8 . When 
the same calculation was performed on the 64-atom cell 
only 20 SCF steps were needed. This favorable scaling 
with respect to both energy cutoff and system size, and 
the inherent stability of multigrid methods offer the pos- 
sibility of performing electronic-structure calculations on 
very large systems. 

In summary, we have developed a methodology for per- 
forming large-scale ab initio electronic-structure calcula- 
tions entirely in real space. The method uses highly- 
efficient multigrid techniques to accelerate convergence 
rates, which is particularly important for ill-conditioned 
systems requiring high energy cutoffs or with large length 
scales. In addition, the multigrid method is readily 
adaptable to parallel computer architectures, and its 
block structure makes it very suitable for O(N) ap- 
proaches. 
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FIG. 1. Variation of the total energy as a carbon atom is 
displaced relative to the grid point along a coordinate axis. 
The grid spacing h is 0.336 bohr. 



TABLE I. Mehrstellen discretization weights in 3D (h is 
the grid spacing). 



Grid point 



6h 2 *A 



6*B 



central 
nearest neighbors 
second neighbors 



-24 
2 
1 



TABLE II. Comparison of perfect crystal and vacancy in 
diamond results obtained using 64-atom supercells [eV]. 



Car-Parrinello a Multigrid Experiment 



Perfect Crystal 

Band gap 

Cohesive energy 
Vacancy 

Formation energy 



4.53 
8.49 



4.53 
8.54 

7.07 



5.50 
7.37 



a 35 Ry pHtoff and Y point sampling. 
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TABLE III. Comparison of GaN results obtained with Ga 
3d electrons in valence. See text [eV]. 





Band Gap 


Cohesive Energy 


Fiorentini et al. a 


2.00 


10.89 


Wright et al. b 


1.89 




Multigrid 


1.79 


10.68 



a Full potential LMTO.E3 p 
b Plane wave calculation with 240 Ry cutoffEJ 
c 32-atom supercell, 160 Ry cutoff and T point sampling. 



TABLE IV. Multigrid convergence tests in an 8-atom su- 
percell: grid spacing [bohr] , equivalent plane- wave cutoff [Ry] , 
total energy [a.u.], and the number of steps to converge the 
density. 



Grid Spacing 


Cutoff 


Total Energy 


SCF steps 


0.421 


25 


-44.87968 


22 


0.336 


35 


-45.03331 


17 


0.280 


60 


-45.06240 


21 


0.210 


110 


-45.06605 


26 



4 



1.0 



0.5 



0.0 



-0.5 



-1.0 



0.0 0.5 1.0 

Displacement along coordinate axis (units of h). 
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